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Abstract 

We show that configuration space techniques can be used to efficiently calculate the complete 
Laurent series e-expansion of sunrise-type diagrams to any loop order in Z?-dimcnsional space- 
time for any external momentum and for arbitrary mass configurations. For negative powers of 
e the results are obtained in analytical form. For positive powers of e including the finite e° 
contribution the result is obtained numerically in terms of low-dimensional integrals. We present 
general features of the calculation and provide exemplary results up to five loop order which are 
compared to available results in the literature. 
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1 Introduction 



The computation of multi-dimensional integrals corresponding to multi-loop Feynman diagrams is a 
necessity for high precision calculation in quantum field theory since perturbation theory remains the 
main tool of most theoretical analysis' within the Standard Model and beyond 0I21E1E]- An important 
ingredient of the algebraic approach to the evaluation of Feynman diagrams is the integration-by-parts 
technique which allows one to derive and analyze recurrence relations for the sets of relevant multi- 
loop integrals [S |. Using the integration-by-parts technique, many diagrams are reduced to simplified 
subsets of master integrals (see e.g. [HI El El El E3 CU C21 HS1 HH ) 

An important subset of master integrals is represented by diagrams of the sunrise topology, the 
so-called sunrise-type diagrams (also called sunset-type, water-melon, basketball, or banana dia- 
grams) |151 1161 I17j . Diagrams of the sunrise topology have been studied quite intensively in the 
past and many of their properties have been known for quite some time |H| I18| IT9 l I13| 120] . Physical 
applications are also numerous as worked out in e.g. |21l 1221 I2U1 1241 125j . This concerns in particular fi- 
nite temperature calculations using effective potentials |26[ I27 | 128 } I29 [ 130] . There are also applications 
in nuclear and solid state physics where methods of quantum field theory are used for the construction 
of low-energy effective theories |31| . The importance of the sunrise topology in physical applications 
alone justifies the strong interest in sunrise-type diagrams jSH ESI EH ESI EH1 EZ1 EH EH1- At the 
same time this topology is a good laboratory for checking the efficiency of new methods of multi-loop 
calculations |4()| I41| H2*|. The relevant master integrals should be calculated as precisely as possible 
within dimensional regularization, the results of which should be expressed in terms of the complete 
Laurent series expansion in the dimensional parameter e = (D — 4)/2 including also positive powers 
of e 43 j. Higher order terms in the e expansion are needed if the sunrise-type diagram is inserted 
into a divergent diagram or when one is using the integration-by-parts recurrence relation which can 
generate inverse powers of e. For example, recently the Laurent series expansion for four-loop master 
integrals has been found using numerical methods [Hj. The numerical evaluation of Feynman dia- 
grams has become a valuable tool in higher order loop calculations because the extreme complexity 
of these calculations often precludes an analytical approach and thus the numerical approach is the 
only possibility to obtain physical results |45| . 

In this paper we develop and describe a configuration space method which allows one to efficiently 
evaluate the coefficients of the Laurent series expansion in the dimensional parameter e for sunrise- 
type diagrams [IB]- The method is simple to use and it works well for arbitrary mass configurations, 
arbitrary values of the external momentum and any number of loops. 



In configuration space, the correlator function described by a sunrise-type diagram with (n + 1) 
lines (corresponding to n loops, cf. Fig. [J) is given by the product of propagators D(x,m), 

n+1 




Figure 1: n-loop diagram with (n + 1) lines of the sunrise-type topology 
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and/or their derivatives if necessary (for details see Refs. |461 1471 The propagator D(x,m) of a 

massive particle with the mass parameter m in D-dimensional (Euclidean) space-time is given by the 
momentum space integral which can be evaluated in terms of Bessel functions, 



, 1 r e - i( P- x )d D p _ (mx) x K x (mx) 
D{x,m) - j—jn J p2 + m 2 ~ ( 2 ^)A+i x 2A <-> 



where we write D = 2A + 2. Here K\{z) is a modified Bessel function of the second kind (see e.g. [49 ). 
In the zero mass case the propagator simplifies to D(x,0) = r(A)/4-7r A+1 x 2A . It is obvious that the 
correlator function II(x) in configuration space contains no integration at all. In this respect it is a 
kind of analogue to a tree diagram in momentum space. A calculation of a sunrise-type diagram is 
necessary only if one wants to calculate the diagram in momentum space. This requires a Fourier 
transformation, 

II n (p) = / U n (x)e i ^d D x. (3) 



Note that the required integrals are basically scalar which makes the angular integration in Eq. (|3"|) 
simple in D-dimensional space-time. One obtains one-dimensional integrals 

U n (p) = 2ir x+1 J™ A J x (px)U n (x)x 2X+1 dx, (4) 

where p = \p\ and x = \x\ are the absolute values of p^ and x^, respecticely, and J\(z) is the Bessel 
function of the first kind. For integrals with additional tensor structure the plane wave function e l ( p ' x ^ 
occurring in the Fourier integral in Eq. © has to be expanded in a series of Gegenbauer polynomials 
Cj(w) ISOlEll- Because this expansion does not change the principal structure of the expressions, the 
representation given by Eq. (0J) is quite universal and will be used in the following. 

The paper is organized as follows. In Sec. 2 we explain the treatment of the singular part of the 
integral while the calculation of the nonsingular part is dealed with in Sec. 3. After a striking but quite 
simple example we present four- and five-loop calculations of the bubble diagram in Sec. 4. Finally, 
in Sec. 5 we give recipes to treat so-called irreducible numerator factors in sunrise diagrams. In Sec. 6 
we present our conclusions. 



2 ^-expansion of the singular part 

The ultraviolet (UV) singularities that appear in the correlator function H(p) are related to the small 
x behaviour of the integrand in Eq. Q. These singularities can be subtracted using the standard 
.R-operation 52 . The simplest way to do this is to expand the relevant Bessel functions at small x 
and then to integrate the resulting integrand in D = 4 — 2e dimensions. However, in order to avoid 
infrared singularities coming from the large x-region of integration, some parts of the unexpanded 
Bessel functions should be retained. From a technical point of view one or two Bessel functions in 
the integrands can be kept unexpanded to provide the necessary cutoff. This is possible because the 
integrals containing one or two Bessel functions in the integrand can be done analytically. Indeed, for 
any //, v one has EI] 

/•CO 

/ x^'~ 1 K v {mx)dx 
Jo 

/ x^~ K v {mx)K v (mx)dx 
Jo 



^rfs+'Vrswsws-*)- ( 5) 
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However, because one or two of the Bessel functions are kept unexpanded, this method breaks the 
natural symmetry between the masses of different lines. An alternative recipe is to introduce a damping 
factor into the integrand with an adequate power of n such as e.g. 

h(x) = e~^ 2 £ = e~^ 2 {l + A 2 + ^ A 4 + A 6 + ■••}• (6) 

One then expands all massive propagators for small x keeping only the singular terms of the integrand. 
For the extraction of the UV poles this is the most convenient way to proceed. At small x the damping 
factor behaves as h(x) = 1 + 0(x 2n ) which does not change the structure of the singularities of the 
integrand at small x if n is sufficiently large (depending on the number of lines and the dimension of 
space-time). Still another possibility is to use a hard cutoff in configuration space by introducing a 
cutoff in the integration at some given xq and to extract the poles from the integral over the finite 
interval (see e.g. Ref. [23). Which of the three methods is the most suitable depends on the problem 
at hand. We found that the damping factor method is most convenient for the calculation of the pole 
parts. In the equal mass case it is advantageous to use the Bessel function method for the calculation 
of the nonsingular parts of the e-expansions since one generates more compact expressions for the 
necessary subsequent numerical evaluation. 

When one performs a series expansion of the integrand near the origin, one easily obtains the 
singular parts for any sunrise-type diagram for any mass configuration. As an example we consider 
the genuine two-loop sunrise diagram with three different loop masses, 

U 2 (x) = D(x, mi)D(x, m 2 )D(x, m 3 ) (7) 

for D = 4 — 2e space-time dimensions (the index "2" in Ii2(x) stands for the two-loop case). In this 
case one needs two UV counterterms for the renormalization of the divergent integral, 

U r 2 cn {x) = n r 2 cg (x) - C 5{x) - C 2 d 2 5(x) (8) 

with d 2 = d^d^ which will lead to a factor —p 2 in momentum space. The counter terms read 

C = M 2 J U 2 (x)d D x = M 2 J D(x,m 1 )D(x,m 2 )D( y x,m 3 )d D x, 

C 2 = Ju 2 {x)x 2 d D x = ^ J D(x,m 1 )D(x,m 2 )D(x,m 3 )x 2 d D x. (9) 

It is clear that the program of UV renormalization requires the calculation of vacuum bubble diagrams. 

Note that we use a special normalization convention for the integration measure. This is the usual 
integration measure for massive vacuum integrals that considerably simplifies the expressions for the 
counterterms (pole parts) by removing some awkward terms containing the Euler constant 7 or tt 2 . 
The normalization factor is given by 

/ lA„\2-e \ n 

(10) 

" c > / 

for the n-loop (or (n + l)-line) sunrise-type diagram. 

Starting with the singular pieces for the genuine two-loop sunrise diagram with equal masses m, 
we obtain 

Go = m 2 {-A-|} + O(e°), C 2 = - £ +0(e°) (11) 

where the renormalization scale [i is set to \x = m. After the counterterms have been determined, 
external lines carrying external momenta can be added because the pole parts of the vacuum bubbles 
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and the sunrise-type diagrams are identical. The generalization of this example to ra-loop sunrise-type 
diagrams with any mass configuration is straightforward. In the main text we list a few sample results 
for equal masses, again setting the renormalization scale /i to fj, = m. This includes sunrise-type 
diagrams with vanishing outer momenta called bubble diagrams. We compute the four- loop bubble 
diagram also calculated by Laporta [H] whose result we reproduce. In the sample results below the 
normalization factor M n is included according to the number of lines/loops that makes the definition 
of Il2(p 2 ) a bit different from that used in Eqs. (|3I4|) . The sample results are 

n!(p 2 ) = - + O(e ), 



£ 



= ™ ! {-i-IK+^°>. 

<W = <0 = ^{| + § + |} + O(£ 0), 

=- 2 ) = ™ 4 {!+§+I}+°(A 

~ , o , fi f 5 35 4565 583451 _ n , 

TTf„ 2 - m 2 ) rr, 6 / 5 45 4255 106147 ] | Q( r ^ 

~ /2 , s f 3 33 1247 180967 8985171 ^. n , 

- /2 2 , 8 f 3 16 49 6967 17060631 n , 

iW = 0) = m 10 ( — _ 133 _ 238 _ 77329 _ 21221921 _ 25963723871 
bV ^ ; I 2e 6 6e 5 3e 4 360e 3 43200e 2 2592000e J ^ ; v ; 

Note the dependence on the external momentum p 2 . This dependence has its origin in the derivatives 
appearing e.g. in Eq. (JSJ). The coefficient of the leading singularity in e is independent of p 2 . In 
Appendix A we list results for unequal mass configurations up to four-loop order. When setting all 
masses equal the results of the general mass case can be seen to agree with the above equal mass 
results. 

After having determined the singular parts of the Laurent series expansion using the damping 
factor method what remains to be done is to calculate the coefficients of the positive powers of e in 
the e-expansion including the finite e° term. In order to determine the nonsingular parts one needs 
to resort to the Bessel function method. What is technically needed is to develop a procedure for the 
e-expansion of Bessel functions. This will be the subject of the next section. 

Starting with the next section we discuss only bubble diagrams which have a simpler structure and 
allow for a comparison with results in the literature (e.g. with Ref. [B]). It is clear that nonsingular 
parts can also be calculated for sunrise-type diagrams with p 2 ^ with our methods. For p 2 > ^ m 2 
the calculation gives rise to absorptive parts represented by the spectral density (cf. Refs. [4*6 l l4*7 | 156] ) 
which will not be discussed any further in this paper. 



3 e-expansion of the nonsingular part 

The nonsingular parts in the e-expansion of sunrise-type diagrams can also easily be calculated with 
the help of configuration space techniques. However, in contrast to the singular coefficients calculated 
in the last section, in the general case the computation of the nonsingular coefficients requires a 
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numerical evaluation. A technical problem which appears in the computation of higher order terms 
of the e-expansion is the necessity to expand the Bessel functions in their indices. For the first order 
in e the relevant corrections are known and can again be re-expressed in terms of Bessel functions. 
Using 



dK u (z) 
dv 



k=0 



for the derivative of the Bessel function K u (z) with respect to its index near integer values of this 
index, we obtain the series expansion 

K- e (x) = K (x)+O(e 2 ), 

K^ E (x) = K 1 (x)--K (x)+O(e 2 ), 

x 

K 2 -e(x) = K 2 (x)--K l {x)-%K Q (x)+0{e 2 ). (14) 

X x z 

for the first few Bessel functions with non- integer indices (for details cf. Ref. [56 ). Note that the 
first two results suffice to find the expansion for the Bessel function with any integer index due to 
the recurrence relations that connect three Bessel functions with consecutive indices. To the best of 
our knowledge one has to proceed numerically for the higher order terms. At least we were not able 
to find a general procedure for the analytical evaluation of the finite parts of the e-expansion. Some 
analytical results can be found in [571 EE] • 

A convenient starting point for the e-expansion is the integral representation 

roo 

K v {z)= e- zcosht cosh(vt)dt (15) 
Jo 

Then the expansion for K_ £ {z) reads 

roo °° F 2n 

K- £ (z)= e- zcosht cosh(-et)^ = £-— f 2n (z) (16) 
Jo ^ (2n)! 

roo 

f k (z) = / t k e- zcosht dt. (17) 



where 



The family of functions fk(z) is rather close to the original set of Bessel functions K v {z) and can easily 
be studied both analytically and numerically. The limits at z — > and z — > oo are known analytically 
and are simple. They allow for an efficient interpolation for intermediate values of the argument. 
For the expansion of the second basic function K\— e {z~) we write 

poo °° F 2n oo 2n+l 

Ki-e(z) = / e~ zc osht co S h(t-et)dt = £ — - -a 2n (z) - £ — — — b 2n+1 (z) (18) 

with 

/■oo roo 

a k {z)= t k e- zcosht coshtdt, b k [z) = t k e~ z cosht sinh t dt. (19) 
jo Jo 

Integration by parts and parameter derivatives can be used to further reduce integrals containing 

any power of the hyperbolic functions sinh(i) and cosh(t) in fk(z). The functions f k (z) allow one to 

calculate the higher order coefficients in the e-expansion. Therefore, the two functions a k {z) and b k (z) 

are again related to the functions f k (z), 

a k {z) = -^fk(z), b k {z) = - z h- X {z). (20) 
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Using these relations, we obtain 



n=0 v ' 

Due to recurrence relations in the index v for the family K v {z) these two formulas suffice to calculate 
the Bessel function K u (z) for any v. 

The functions fk(z) satisfy the differential equation 

( d 2 1 d \ , x fc(fc - 1) 

(a? + ;s " ') A(2) - -V^w < 22 > 

related to the Bessel differential equation. The small z behaviour 

AW-STT ta * H G)( 1 + (s5j)) (23) 

can be found by using yet another representation for the function fk(z), obtained from Eq. ()17j) by 
the substitution u = cosht, 



f k {z) = / -== ln*(u + Vtf=T)du (24) 
or directly from the behaviour of the function K u (x) at small x, 

K- e (x) = -sinh(eln(l/^)) . (25) 

In Fig-Elthe functions f n ( z ) ar e compared with Kq{z) for various values of n. Note that fo(z) = Kq(z) 
and Oq(z) = —Jq(z) = K\(z). All curves are very smooth and their functional behaviour is similar to 
that of the original Bessel function. 

While in four-dimensional space-time the massive propagator contains the Bessel function Ki_ e (z), 
for D = 2 dimensional space-time the basic function is the Bessel function K— e {z), since the propagator 
reads 

D(x,m)\ D=2 = ^-K (mx). (26) 

As an example for the numerical calculation of the e-expansion using Bessel functions we consider a 
toy model integral related to the one-loop case in two dimensions. We select this example because the 
integral is finite and analytically known, so that we can compare our numerical calculation with the 
exact answer. Using Eq. ©, we obtain 

2 J K_ £ {x)K^ £ {x)x dx = T(1 + E )r(l - e). (27) 

The expansion in e is given by 

r(i - e )r(i + e) = = i + + + ( e 6) (28) 

sin(7T£) 6 360 

On the other hand, we can use the expansion 

2 4 

K^{x) = f (x) + £ -f 2 (x) + £ -U{x) (29) 
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Figure 2: Comparison of functions f n (z) for different values of n, plotted on a logarithmic scale 



to rewrite the integral in the form 



2 K- £ (x)K_ £ (x)xdx = 2/ f (x) 2 xdx + 2e 2 f (x)f 2 (x)xdx + 
Jo Jo Jo 

-4 r oo 



pi roc pt roo 

fo(x)h(x)xdx + - f 2 (x) 2 xdx + 0(s i 
o Jo 2 Jo 



(30) 



Using the explicit expressions for the functions fk we checked by numerical integration that the 
identities 



roo 

2 / Uxfxdx = 1, 
Jo 
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2 / f (x)f 2 (x)xdx 
Jo 

(fo(x)h(x)+3f 2 (x) 2 )xdx 



TV 

~6~' 

7tt_ 
360 



(31) 



are valid numerically with very high degree of accuracy. We have implemented our algorithm for 
the e-expansion of sunrise-type diagrams as a simple code in Wolfram's MATHEMATICA system for 
symbolic manipulations and checked its work-ability and efficiency. 

4 Examples 

In this section we compute some further examples using our techniques. First we check on a recent 
result obtained by Laporta for sunrise-type four-loop bubbles [44] . In Ref. [11] the difference equation 
method is used to numerically obtain the coefficients of the Laurent series expansion of all four-loop 
bubble master integrals. Our results for the sunrise-type topology derived by using configuration 
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space techniques provide an independent check for the results in Ref. |44| where momentum space 
techniques were used to calculate the whole set of four-loop bubble master integrals. This check may 
help in establishing further confidence in the results of Laporta. As a new application we present the 
five-loop result for the scalar bubble topology. 

4.1 Four-loop vacuum bubble 

The four-loop quantity of interest (Vi in the notation of 01]) is the master integral 

n 4 (p 2 = 0) = Ma J D(x, mfd D x (32) 

where the normalization factor A/4 is defined in Eq. (jlOjl . The Laurent series expansion of the four-loop 
master integral has been calculated numerically by the difference equation method with high precision 
in [II]. In Sec. 2 we explained how to obtain the coefficients of its singular part. In fact, we wrote 
down explicit results for the singular part of the four-loop master integral. In this section we shall 
compute the finite part and the first three coefficient functions of its e-expansion numerically using 
configuration space techniques. 

The general idea is really quite straightforward. As explained before, UV divergences reveal 
themselves in configuration space as singularities of the integrand at small x. We subtract these 
singularities, obtain a regular integrand, expand the integrand in e and then finally integrate the 
integrand numerically 47 . 

Let us describe our procedure in more detail. We write 

n 4 (0) = A/4 J D(x, mfd D x = A/4 J D(x, m){D(x, m) - A(x) + A(x)fd D x 
= A/4 J D(x, m) ((£>0, m) - A(x)) 4 + 4(D(x, m) - A(x)) 3 A(x) 

+6 (D(x, m) - A(x)) 2 A(xf + 4 (D(x, m) - A{x)) A{xf + A(j;) 4 ) d D x (33) 

with 

A(x) = D(x,0) + £>i(:r,0) + D 2 (x,0) (34) 

where the functions D(x,0), Di(x,0), D2(x,0) subtract singularities of D(x,m) at small x. Their 
explicit expressions are given in Appendix B. In fact, these functions represent a formal expansion 
of the massive propagator D(x,m) in the mass parameter m for small m since the real expansion 
parameter is the dimensionless quantity mx. 

One propagator factor D(x, m) in the integrand is left unsubtracted. At large x it provides an IR 
cutoff. The last two terms of the integrand in Eq. ()33|) can be integrated analytically. They contain all 
poles in e and are expressible through Euler's T-functions. As expected, the pole part coincides with 
the expression in Eq. (|T2*|). Because the analytical expression (as given up to order e 3 in Appendix C) 
is rather lengthy, we present only its numerical evaluation, 

n^ na (0) = m 6 ( - 2.5e~ 4 - 11.6666667e~ 3 - 31.701389e~ 2 - 67.5289356" 1 

-15871.965743 - 142923.10240e - 701868.64762c 2 - 2486982.5547e 3 + 0(e 4 )) . (35) 

The first three terms in Eq. ()33|) can be integrated numerically for D = 4 (i.e. no regularization is 
necessary) since they are regular at small x. The analytical expression for the functions to be integrated 
is rather long. The zeroth order e-coefficient is found in Appendix D (see also the discussion of the 
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Figure 3: Integrands for numerical integration in case of the four-loop bubble at different orders in e 
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integration procedure given in Appendix D). However, as shown in Fig. |2J the functions themselves 
show a very smooth behaviour which renders the numerical integration rather simple. We obtain 

f[™ m (0) = m 6 (l5731. 745122 + 142349.56687e + 699112.42072e 2 + 2468742.6339s 3 + 0(e 4 )) . (36) 

The sum of both parts gives 

n 4 (0) = m 6 ( - 2.5e" 4 - 11.6666667e- 3 - 31.701389e" 2 - 67.528935e" 1 



-140.220621 - 573.53553e - 2756.22690e 2 - 18239.9208e 3 + 0(e 4 



(37) 



which confirms the known result jll] . The difference to the results of |44| is within the accuracy of the 
numerical evaluation of the integrals. Therefore, we provide an independent check of this important 
quantity in full agreement with the results of Laporta. 

It is not difficult to extend the analysis to higher orders in e or to a larger number of significant 
digits in the coefficients of the e-expansion. However, since the technique is rather straightforward 
and simple we do not consider it worthwhile to extend the calculations into these directions. If the 
need arises, the potential user can tailor and optimize his or her programming code to obtain any 
desired accuracy and/or order of the e-expansion. In our evaluation we have used standard tools 
provided by the MATHEMATICA package which allows reliably control the accuracy of numerical 
calculations. Even at this early step of improvement it is obvious that our algorithm is extremely 
simple and efficient. 

4.2 Five-loop vacuum bubble 

In this subsection we present results for the next order p 2 = sunrise-type diagram, the five-loop 
vacuum bubble. We have chosen to extend our calculation to the five-loop case since there exist no 
results on the five-loop bubble in the literature. The integral representation of the five-loop bubble is 
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given by 

fl 5 (p 2 = 0) = AT 5 J D(x,mfd D x. (38) 

Evaluating numerically the analytical part one obtains 

ri^ na (0) = m 8 (3e~ 5 + 16.5e~ 4 + 51.95833e~ 3 + 125.6715e" 2 + 259.98766" 1 

-1360392.5934 - 16888723.177e - 111392297.46e 2 - 518606741. le 3 + 0(e 4 )) (39) 

while the numerical integration of the nonsingular part gives 

fig um (0) = m 8 (l360739.9485 + 16886269.683e 

+111360751.91e 2 + 518295438.0e 3 + 0(e 4 )) . (40) 

The sum of both contributions is given by 

fi 5 (0) = m 8 (3e- 5 + 16.5e" 4 + 51.95833e" 3 + 125.6715e~ 2 + 259.9876e" 1 

+347.3551 - 2453.494e -31545.55e 2 - 311303.1e 3 + 0(e 4 )) . (41) 

One observes huge cancellation effects between the terms obtained by the analytical calculation and 
the numerical integration. Apparently the subtraction procedure chosen here is non-optimal. As men- 
tioned before, the subtraction procedure should really be optimized for any given problem in order to 
avoid a necessity to retain high numerical precision at intermediate steps of calculations. Nevertheless, 
our nonoptimized simple subtraction procedure already works quite reliably with available standard 
computational tools. 

In this section we have described how the configuration space technique works for the case of a 
trivial numerator. However, our method is applicable and quite efficient also for nontrivial numerator 
factors as shown in the next section. 

5 Irreducible numerator factors in sunrise-type diagrams 

For vacuum bubbles of a given topology there can be more than one master integral. For example, 
in the case of the four-loop sunrise-type bubble topology Laporta identified a second master integral 
V2 which has a nontrivial numerator factor which cannot be further reduced [11]. In the momentum 
representation the non-trivial numerator factors contain scalar products of loop momenta which cannot 
be further canceled against denominator factors. To be precise, in the case of a n-loop vacuum bubble 
with (n + 1) massive lines, the numerator factor is trivially reducible for n < 3. However, starting 
with n = 3 a numerator factor is no longer directly reducible in the general case. 



5.1 Three-loop vacuum bubble with irreducible numerator 

As an example let us consider the case of a numerator (ki ■ kj) for a vacuum bubble with four massive 
lines. In momentum space the integral of interest reads 

fWru - f^ D / 2 rd- n/9\Y 3 f (fci • k 3 )d D k 1 d D k 2 d D k 3 

An expansion of the numerator in terms of invariants, 

2(h ■ k 3 ) = k\ + k\ - (h - k 3 ) 2 (43) 
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contains a structure (k\ — k 3 f which is absent in the denominator. Therefore, the numerator is not 
directly (trivially) reducible to unity. 

In configuration space it is not difficult to treat such a non-trivial numerator since one can reduce 
the numerator to derivatives of the propagators. The main features of the reduction method are more 
easily discussed in terms of the equal mass case. The generalization to nonequal mass case is evident. 
In the above example one can write (using m\ = m 2 = . . . = m) 



fj*(0) = 7V 3 J D(x, mf (df.Dix, m)) (d^D^x, m)) d D x (44) 

The integration-by-parts identity takes the form 

J (D(x, m) • • ■ D(x, to)) d D x = (45) 

which, as usual, should be used with the necessary caution. 4 A further useful identity is given by 

J d 2 (D(x, m) • • • D(x, m)) d D x = 0. (46) 

The identities in Eq. (|45|1 and (|46|) lead to a relation between integrals with derivatives. In our example 
involving four propagators one obtains 

J D(x, mf {d^D(x, m) (d»D(x, m))) d D x = -^J D(x, m) 3 d 2 D(x, m)d D x. (47) 

For the last integral we use (— d 2 + m 2 )D(x,m) = 5(x) to end up with 

D( X , m ?e>D(x, m)i ° x = m ' / D (x, ra ) V>» - D (0, m f. (48) 

The value of D(0,m) can be found by integration in momentum space, 5 



d D p 2tt d / 2 f°° p°- l dp 



m 



D-2 



D(0, m) - (2 ^ )D J m2+p2 - (27r) D r(jD/2) J o p 2 + m 2- {47t) d/2 r (! D / 2 )- ( 49 ) 
Finally we have 

ns(0) = -^fia(O) + \NzD{Q, mf (50) 

where 113(0) is a three- loop sunrise- type bubble diagram without numerators, as calculated earlier. 
This relation can explicitly be checked. 

5.2 Four-loop vacuum bubble with irreducible numerator 

As a more realistic example we consider a four-loop diagram that appears as a second independent 
master integral V2 of the sunrise topology in the classification of Ref. |44j . In momentum space the 
second master integral V% has the additional numerator factor {k\ ■ k£f as compared to the first master 
integral V\ with unity in the numerator. The second master integral reads 

n*(0) = U D/2 T(3-D/2)Y 4 x (51) 



(fci • k 4 fd D k 1 d D k 2 d D k 3 d D k 



4 



K + fc?)(mj + (fc2 - k^iml + (k 3 - k 2 ) 2 ){m\ + {k 4 - fc 3 ) 2 )(m| + k\) ' 



4 Why this caution is necessary is illustrated by the simple example of a massless propagator. The massless propagator 
satisfies the equation — d 2 D(x,0) = S(x) where d 2 = 9 M <9 M . Integrating the equation over the whole space and dropping 
the total derivative on the left hand side one obtains the contradiction = 1 (for details cf. |59)'l. 

5 D(0, m) can also be obtained by taking the limit x — > on the right hand side of Eq. @. 
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Turning again to the equal mass case and using the configuration space representation, this integral 
can be written as 

ft*(0) = A/ 4 J D{x,m) 3 (d^d„D{x,m)){d^d l ,D{x,m))d D x. (52) 

It is apparent that by using integration-by-parts techniques this integral cannot be reduced to scalar 
integrals and/or integrals containing d'Alembertians. The easiest way to evaluate such an integral is 
to compute the derivatives directly. This is done by using 

\lz ( z ~" K "W = ~ {*- v - X K V +x(*j) ■ (53) 

This relation can be iterated and gives results for arbitrary high order derivatives of Bessel functions 
K\(z) in terms of the same class of Bessel functions with shifted indices and powers in z. For the first 
derivative we obtain 

m 2A+2 K\ + i(mx) 
; (2vr) A+1 (mi) 

Since the resulting analytical expression for the given line of the diagram lies in the same class as 
the original line, the procedure of evaluation of the integral is similar to the usual one. However, we 
cannot use the second derivative 



8 D(x m) - -x — " A+1V "^ f54) 



m 2X+2 ( . . x^x u 



d^d u D(x,m) = (27r) A+i (mx) A+i [V%H " ^K x+1 (mx) j (55) 

directly under the integration sign. The reason is that if the propagator is considered as a distribution, 
there is a (5-function singularity at the origin which is not taken into account in Eq. (|55|) . Indeed, 
contracting the indices /i an v in Eq. (|55|) one obtains 

d /1 d fl D(x,m) = m 2 D(x,m) (56) 

while the correct equation for the propagator reads (— d 2 + m 2 )D(x,m) = 5(x). Thus, the straight- 
forward evaluation of derivatives is always valid only for x = 0. The behaviour at the origin (x = 0) 
requires a special consideration. In practice, to treat this situation one should not use higher order 
derivatives but stay at the level of first derivatives. 

In order to deal with this situation we introduce another master integral 

f[**(0) = A/" 4 J D(x, m)d^D(x, m)d v D{x, m)d^D(x, m)d u D(x, m)d D x. (57) 

The relation between the two master integrals II|(0) and II|*(0) is found to be 

n^(0) = 3nT(0) - Jm 4 n 4 (0) - lm 2 M^D(0,m) A . (58) 



The quantity II|*(0) can be calculated with the explicit use of first order derivatives within our 
technique. We obtain the analytical result for the pole part 

Il**(0)-m w ( A JZL E^L 4936643 A 

n 4 (0) - m ^- g£4 - i44e3 - mi2£2 - 4i472Qe + 0[e )j (59) 

and the e-expansion in the form 

ri**(0) = ml °( ~ 0.375e -4 - 1.923611e~ 3 - 5.474103e~ 2 - H.90356^ 1 

-27.99303 - 104.5384e - 663.6123e 2 - 3703.241e 3 + 0(e 4 )) . (60) 
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Since the results for 114(0), f[|*(0), and -D(0, m) 4 are known, Eq. (|58|) can be used to obtain the final 
result for the original integral, 



which again verifies the result given in Ref. |44j . 

Differentiation of the massive propagator leads to expressions of a similar functional form which 
makes the configuration space technique a universal tool for calculating any master integral of the 
sunrise topology. This technique is also useful for finding master integrals. Indeed, new master 
integrals appear when there is a possibility to add new derivatives into integrands which cannot be 
eventually removed by using the equations of motion or integration- by-parts recurrence relations. But 
once again: without explicit inclusion of the <5-function only one derivative is allowed. Otherwise one 
runs into problems not seeing some parts (tadpoles) of the result. Therefore, the new master integral 
should contain just one derivative for each line excepting one line. For instance, in the five-loop case 
(six propagators) there will be only two master integrals. 

6 Summary and conclusions 

We have suggested a new efficient technique to compute diagrams of the sunrise-type topology with 
any number of loops at any order of the e-expansion for any mass configuration. For a given number 
of loops the sunset-type topology constitutes only a small fraction of the whole set of multi-loop 
topologies that need to be calculated. Nevertheless, our results on the subset of sunrise-type diagrams 
provide a necessary check on multi-loop results calculated by other techniques and therefore can be 
very useful for many multi-loop calculations. We have also worked out a few examples with nontrivial 
numerator factors. Valuable by itself, our method can be used to check the results of other techniques. 
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A Singular contributions for arbitrary masses 

In the following we present complete results for the singular parts of sunrise-type diagrams with 
arbitrary masses up to four-loop order. The results are given in the MS-scheme in the Euclidean 
domain. 



n2(o) 



m 




(61) 




1 




1 

2? 
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n|(p,mi,m2,m3,m 4 ) = ^3 m 2 m 2 + I p 2 ^ w 2 - ^ m 4 + 2 ^ mfm 2 (4 - 3^) ] + 



Ut(p, mi ,m 2 , m3 ,m 4 ,m 5 ) = -JL £ m 2 m 2 m 2 + 

p 2 mfm 2 — y^^(mfm 2 + mfm 4 ) + 2 ^ m 2 m 2 m\(h — 4^) + 



48e 3 
1 

288e 2 



^ (2p 4 2 m 2 - 6p 2 2 m 4 + 2 £ m? + 3p 2 ]T m 2 m 2 (ll - 8£) + 
-6^(m 4 m 2 + m 2 m 4 )(7-4^) + 12 ^ m 2 m 2 m|(15 - 20l { + 2^ 2 + 6^) j + 

"172^ ( 3p6 + ? ^ (35 " 24 ^ } " 18p2 ? ^ (21 " 8 ^ + 2 ? ™* (7T " 2Ui) + 

\ i i i 

+9p 2 m 2 m 2 {7l - 88^ + 8^ 2 + 24^) - 216 ^(m 4 m 2 - m 2 mf)ti + 
- 18 ^(m 4 m 2 + m 2 mp (49 - 56£ + 4^ 2 + 12^) + 

+24 ^ m 2 m 2 m 2 k (W5 - 180^ + 30^ 2 + 90^j - 2^ 3 - 18^- - 12^4) J (62) 

where £i = m(m 2 //U 2 ). The indices i, j, and k run over all mass indices. One can check that the 
general results listed in this Appendix reproduce the equal mass results listed in the main text. 

B Explicit form of subtraction terms for the small x singularities 

The leading singularity at small x is given by the massless propagator of the form 

^'°) = 4^A- ( 63 ) 
The next order of the small x-expansion for the propagator D(x,m) is explicitly given by 

^)=^(f§m-fsfffi^v (64) 



Att x+1 x 2X \\2j 1-A \2J A 

This term is suppressed relative to the first term by one power of x 2 at small x in four-dimensional 
space-time (however, this is not the case for two-dimensional space-time with A = 0). The term 

*m> - ^ (i y (d) 2 . (.) - e|^i) m 

is further suppressed by one power of x 2 at small x. Therefore, the full subtraction of the three terms 
gives a rather smooth behaviour at small x which is sufficient to obtain a regular integrand for the 
numerical integration. 
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C Analytical results for the four-loop sunrise diagram 



In this appendix we present some more details of our calculations for the four-loop sunrise diagram. 
For the analytical evaluation we take the last two terms of the integrand from Eq. ()33j) . One has 
to integrate a product of two Bessel functions with powers of x which can be done analytically. The 
explicit expression for the e-expansion of that part of the integral which is evaluated analytically reads 



nf a (o) = m 6 - 



5 35 4565 58345 



2e 4 3e 3 144e 2 864e 



1456940638037 17099tt 2 3857tt 4 2525968C(3) 

7779240000 24 10 + 105 

/ 55171475321621447 2457509vr 2 1292537vr 4 
+ \ 1633640400000 + 144 175 

+ 6752^C(3) +1653fa2c(3)+59508c(5) j e 

/ 10610679621089130529 92781949vr 2 4290113759vr 4 22591vr 6 
+ \ 68612896800000 + 864 110250 14 

+ 9524 ™ c(3) + 244476 ' 2 « 3 > - + ^^y 

( 5963907632629558995931 1325204033tt 2 
+ \ 14408708328000000 + 5184 

464379085699vr 4 48529231vr 6 



6615000 2205 
312138383154103C(3) 7285043tt 2 C(3) 43529 4^ _ 238229084C(3) 2 
277830000 + 6 + % 105 

+ 135830 2 1 9 1 4 2 Q 97C(5) + 247950vr 2 C(5) + 1190160C(7)) e 3 + O(e) 4 } . (66) 

This expression shows the real complexity of the calculation which reveals itself in the structure of the 
results. The main feature is that the terms cannot be simultaneously simplified to all orders in e. By 
a special choice of the normalization factor one can make the leading term and, in fact, even all pole 
terms simple, but then the higher order terms contain rather lengthy combinations of transcendental 
numbers that are not reducible in terms of standard quantities such as the Riemann ^-functions. Note 
also that the rational coefficients of transcendental numbers are very big and there is a huge numerical 
cancellation between the rational and transcendental parts of the answer (see also the discussion in 
Ref. 56 ). 



D Integrand for numerical integration 

For the numerical evaluation we take the first three terms of the integrand from Eq. ()33|) . One has 
to integrate them numerically as there is a product of three or more Bessel functions which is too 
complicated to be done analytically. To find the e-expansion of the integral one has to first expand 
the integrand in e. The expression for the e-expansion is quite lengthy. We therefore give explicit 
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results only for e = 0. For this part the integrand for the numerical integration over z = mx reads 



ti( , ,.,2 384 384 768/ 24 480/ 576/ 2 
l (x) = m 6 66 - 108/ - 144/ 2 + 192/ 3 + —= r + — j- + =- + 



z 6 z 4 z 4 z 2 z 2 z 2 



+ 117/ V + 24/ 3 z 2 + 24/V + 



16 2 32 

201/.2 4 9/V „o , l4 4 75z 6 405/z 6 531/ 2 z 6 

H 1 24/ 3 z 4 + 12/ 4 z 4 H h 

16 2 512 128 

3J3 Q 7 K,8 oor;^ qic:/2^8 M/3,,8 



15/ 3 z 6 9/ 4 z 6 375z 8 825/z 8 315/ 2 z 8 51Z 3 

+ — T- + "7777777 777777- + 



2 4 2048 1024 256 64 



18752 10 375/z 10 225/V 15/ 3 z 10 3/V° 

+ TT77777777- 7777777" + 777777 77777" + 



131072 8192 4096 512 512 

-512 384 768/ 24 288/ 384/ 2 
— + 

z z z 




-m° I — =— + — f - + — + - — — - 52z + 1201 z - 64/ 3 z 



— - 21/z 3 + 48/V - 24/V + — - + 9/V - 3/V 

8 32 16 

125z 7 75/z 7 15/ 2 z 7 iV\„,. 2 fi 128Ki(z) 5 

+ + \Ki{z) 2 + m & 67 

512 128 32 8 I u ' z 2 



Here / = ln(e 7 z/2), z = mx, and 7 = — r'(l) is Euler's constant. As shown in Fig. |S1 the plot of this 
function as well as the shapes of the corresponding functions in higher orders of e are very smooth 
and quite similar. The analytical expressions for higher orders in e, however, become much longer. 
Note that the new functions f n {z) first appear at order e 2 . 

The smoothness of the zeroth order integrand as shown in Eq. 1)67(1 implies that the numerical 
integration is quite easy to execute. Because the integrand vanishes exponentially for large values of 
z and has no singularities of the kind z In z for small values of z, the integration can in principle range 
from to 00. However, for practical reasons we had to instruct MATHEMATICA (which we used 
for all of our calculations presented here) that the integrand vanishes for z = 0. On the other hand, 
the asymptotic expansion of the integrand together with the integration measure is dominated by the 
term 

^fz 10 lnV*/2)tfi (^ 3 (Kx(z) - ^e- z for z -+ 00). (68) 
Integrated over z from A to 00, this part gives a contribution 

^^A 25 / 2 lnVA/2)e- A (69) 
512 

and terms which are of subleading order. Therefore, the result can be well estimated by 

nf m (0) = 2tt 2 / UT m (x)x 3 dx « 2vr 2 / UT m (x)x 3 dx + ^J!La 25 /2 i n 4 ( e TA/2) e - A (70) 
Jo Jo 512 

and A can be adjusted in such a way that any desired precision is obtained. 

A possibility to avoid any kind of cutoff is to change the integration variable in the sense that the 
interval [0, 00] is mapped onto [0, 1]. Then the integration can be done numerically with the additional 
information that the integrand vanishes identically at both end points. Possible transformations of 
this kind are for instance z = ln(l/t) or z = (1 — t)/t for t E [0, 1]. 
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